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Recent theoretical and experimental investigations have demonstrated the role of certain invariant manifolds, 
termed burning invariant manifolds (BIMs), as one-way dynamical barriers to reaction fronts propagating 
within a flowing fluid. These barriers form one-dimensional curves in a two-dimensional fluid flow. In prior 
studies, the fluid velocity field was required to be either time-independent or time-periodic. In the present 
study, we develop an approach to identify prominent one-way barriers based only on fluid velocity data 
over a finite time interval, which may have arbitrary time-dependence. We call such a barrier a burning 
Lagrangian coherent structure (bLCS) in analogy to Lagrangian coherent structures (LCSs) commonly used 
in passive advection. Our approach is based on the variational formulation of LCSs using curves of stationary 
“Lagrangian shear”, introduced by Farazmand, Blazevski, and Haller [Physica D 278-279, 44 (2014)] in the 
context of passive advection. We numerically validate our technique by demonstrating that the bLCS closely 
tracks the BIM for a time-independent, double-vortex channel flow with an opposing “wind”. 


Looking into a pocket near the edge of a lazy 
stream, one may see froth and small plant partic¬ 
ulates twisting and evolving in sinuous strands. 
While streams are always subject to change—a 
branch falls in, a rock rolls, a fish swims nearby— 
the observed pattern changes little due to these 
differences in the stream’s current. An explana¬ 
tion for this is that the pattern organizes around 
a robust set of curves that form a “skeleton” for 
particulate advection. 

While it is clear to the eye that such curves ex¬ 
ist, it is another thing to define them precisely. 
A variety of techniques have been developed for 
this purpose, and while the original techniques 
were limited to flows that were either constant 
or periodic in time, the most general modern 
techniques can analyze flows with arbitrary time- 
dependence. 

Considering now not a stream, but the entire 
ocean, one can find a teaming algae population 
within large eddies. The “living pattern” created 
by the swirling collection of replicating organisms 
is different from the passive particulates because 
the algae cluster is active—replication increases 
the size of the algae patch. There is a corre¬ 
sponding new set of “active” curves underlying 
the pattern formed by this growing front. 

In this paper, we develop a technique to extract 
the most attracting and repelling curves for such 
active material within an arbitrary flow defined 
over a finite time interval. We imagine eventual 
application to improving management of oceanic 
algae blooms or to informing the design of com¬ 
bustion systems. 
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I. INTRODUCTION 

A growing interest has recently developed in fluid sys¬ 
tems that combine the phenomenon of (chaotic) advec¬ 
tion with that of front propagation. This interest is 
due in part to geophysical and technological applica¬ 
tions, such as the growth of plankton blooms in oceanic 
flowflM], the growth of the ozone hole in the atmo¬ 
sphere, and the spreading of chemical reactions in mi¬ 
crofluidic devices. Another important driver of this field 
has been a steady accumulation of experimental ob¬ 
servations and measurements of (Belousov-Zhabotinsky) 
chemical reaction fronts spreading in driven, laboratory- 
scale flowiPHj^, as well as early theoretical predictions on 
such advection-reaction-diffusion systemsEKS. 

A recent advance in understanding advection-reaction- 
diffusion (ARD) systems is the realization that invariant 
manifold theory, familiar in the analysis of passive advec- 
tive mixingEDQiD is also applicable to the analysis of front 
propagation in fluids^® 25 '. Such manifolds have been 
called burning invariant manifolds , or BIMs, to distin¬ 
guish them from their analogous, but distinct, advective 
counterparts. The most important property of BIMs is 
that they form time-invariant, or time-periodic, one-way 
barriers to front propagation in fluid flows that are them¬ 
selves either time-invariant or time-periodic, respectively. 
The unstable BIMs are attracting structures, in that an 
initial stimulation near such a BIM ultimately creates a 
reaction front in the fluid that converges upon the BIM. 
In the most striking cases, this front persists for arbitrar¬ 
ily long times, forming a frozen or pinned front, whose 
profile follows the BIM^S. Alternatively, stable BIMs are 
repelling structures. Initial stimulation points on either 
side have qualitatively distinct future behavior. Various 
other aspects of passive invariant manifolds extend to 
BIMs. For example, BIMs lead to a theory of turnstiles 
for front propagation in time-periodic fluid flowJ? 

Despite the successful application of BIMs to ARD sys¬ 
tems, their applicability has been limited to flows that are 
either time-independent or time-periodic. The present 
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paper addresses this limitation by developing a theory 
of prominent, one-way barriers to front propagation in 
general fluid flows specified over a finite interval of time, 
with no assumption on the time-dependence over this in¬ 
terval. These barriers, called burning Lagrangian coher¬ 
ent structures (bLCSs), are either locally most repelling 
or most attracting. We validate this theory by numeri¬ 
cally demonstrating that the bLCS reproduces the BIMs 
for time-independent flows considered over a sufficiently 
long time interval. A future publication will assess this 
theory in the general case of unsteady flows. 

Our work is motivated by recent progress in the use 
of finite-time Lyapunov exponent (FTLE) fields^®!) anc i 
the associated Lagrangian coherent structures (LCSs) 
for passive advection in unsteady flowflSHMl A com¬ 
mon approach in advection studies is to assume that 
ridges of the (forward) FTLE field are (repelling) LCSs 
that approximately separate trajectories with different 
future behavior. The FTLE ridge approach has its defi¬ 
ciencies, however, and can produce both false negatives 
and false positives when identifying LCSsP^, so it is best 
used, perhaps, as an intuitive diagnostic tool. That said, 
if one were to implement the FTLE-ridge approach for 
front-propagation dynamics, one would immediately en¬ 
counter a challenge of dimensions. This is because the 
front-element dynamics in two spatial dimensions is nat¬ 
urally represented as a three-dimensional dynamical sys¬ 
tem, with x and y the two spatial coordinates and 8 the 
orientation of the front element. (See Sect. |TTJ) Note 
that the objects of our interest, the bLCSs, are still one¬ 
dimensional curves, just like the BIMs they seek to gener¬ 
alize; that is, a bLCS is not a codimension-one object in 
our study. Thus, a possible implementation of the FTLE- 
ridge approach to the front-element dynamics would not 
seek to follow ID ridges in a 2D space, or even 2D ridges 
in a 3D space, but rather to follow ID ridges in a 3D 
space. We see no obvious reason why such curves in the 
3D phase space would correspond to fronts in the cry- 
position space, since a curve in the 3D space must satisfy 
the front-compatibility criterion (see Sect. 0 for it to 
be the lift of a front in cry-space. It remains an open 
question, then, whether the ID FTLE ridges would even 
yield fronts in the cry-position space, and assuming they 
did, these would undoubtedly suffer the same deficiencies 
noted in Refill for passive advection. 

Instead of seeking FTLE ridges, we base our approach 
on the recent work of Farazmand, Blazevski, and Haller 
(FBH) 36 !, which uses a variational principle to define the 
LCS. The advantage of a variational formulation is that it 
is based from the start on a search for curves in cry-space. 
These curves could be material lines, as in FBlPI or they 
could be fronts, as in the current study. The fact that the 
evolution of each front element depends on its orienta¬ 
tion 8 is naturally accommodated within the variational 
approach. Following FBH^I, we base the variational in¬ 
tegral on the Lagrangian shear, generalized to the front 
propagation dynamics. The output of the variational ap¬ 
proach, when combined with the front compatibility cri- 



FIG. 1. (color online) Each front element evolves according 
to Eqs. Q. 


terion, is a flow on a 2D surface, which we call the shear¬ 
less surface, embedded within the 3D phase space. This 
flow foliates the shearless surface into shearless fronts, 
which are candidates for a bLCS. This reduction to a 
2D surface tames the challenge of dimensions discussed 
above. One can, for example, plot burning FTLE (bF- 
TLE) fields on the 2D surface. However, the dimensional 
challenges have not been completely eliminated, since the 
shearless surface can have a complicated geometry within 
the 3D phase space, with folds and intersections that are 
not easily represented in a 2D projection. This geometry 
raises new considerations in the search for bLCSs that do 
not arise when searching for LCSs in passive advection. 

This article has the following structure. Section |TT] 
reviews our prior dynamical systems approach to front 
propagation in fluid flows and introduces the concept of 
BIMs as one-way barriers. Section m contains the core 
derivations. The main result is Theorem[2]of Sect. |III G[ 
which identifies perfect shearless fronts as the candidates 
for a bLCS. Sect. |III H| describes how to compute these 
shear less fronts, and Sect. |III T| formulates the selection of 
the (repelling) bLCS from the set of shearless fronts as 
finding that front which maximizes the average normal 
repulsion. Section [TV] applies the bLCS approach to a 
time-independent double-vortex in a channel with wind. 
We study both short (Sect. IV B) and long (Sect. IV C) 
evolution times to show that the bLCS closely reproduces 
the BIM of the time-independent flow. 


II. PRELIMINARIES 

Following previous studietP^II, we give a brief account 
of front propagation from the dynamical systems per¬ 
spective. We make the sharp front assumption, in which 
the front marks a clear delineation between a reacted 
or “burned” region and an “unburned” region. Then, a 
front in a two-dimensional fluid flow can be represented 
as a curve R(s ) = (r(s),0(s)) in the three-dimensional 
xyd- space, where s is an arbitrary smooth parameter¬ 
ization. Here 8 denotes the orientation of the normal 
n = (sin0, — cos0) and tangent g = (cos 8, sin 8) to the 
front, where n x g = +1. (Fig. 0 ) The normal vector 
points in the propagation, or burning, direction of the 
front. 

Not all curves in xyd-sp&ce represent fronts. There 
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must be pointwise agreement between the tangent vector 
r' = dr/ds of the curve r(s) and the orientation g(8) 
given by 0. We call this the front-compatibility criterion 

rj 7 = ±g = ±(cos#,sin#). (1) 

The front propagation speed vq in the local fluid 
frame, which we call the burning speed , is assumed to 
be isotropic and homogeneous throughout the fluid. Fur¬ 
thermore, we assume that Vq is independent of the local 
curvature of the front or any other front property, i.e. vq 
is entirely constant in our analysis. Then each front ele¬ 
ment (r, 0) evolves independently of its neighboring front 
elements according to (see Fig. [l]) 


dr 

a = u + ''" n - 

(2a) 

de 

(2b) 


where u = u(r, t) is the fluid velocity field, dot denotes 
the time derivative, comma denotes differentiation with 
respect to r,, and repeated indices are summed. The 
fixed points of this dynamical system are called burn¬ 
ing fixed points. For time-independent or time-periodic 
fluid velocity fields u(r, t), burning fixed points may have 
any combination of stable (S) and unstable (U) direc¬ 
tions: SSS, SSU, SUU, UUU. The one-dimensional stable 
and unstable manifolds of SUU and SSU burning fixed 
points, respectively, play a special role and are called 
burning invariant manifolds (BIMs). BIMs satisfy the 
front-compatibility criterion and may thus be viewed as 
(virtual) fronts within the fluid. As such, they form one¬ 
way barriers to front propagation, owing to the fact that 
no front can catch up to and surpass another front mov¬ 
ing in the same direction (what we call the no-passing 
lemma 2 A) BIMs are thus natural, dynamically defined, 
geometric objects that restrict the propagation of fronts 
in an advecting fluid flow that is either time-independent 
or time-periodic. 

III. COHERENT STRUCTURES AS SOLUTIONS TO 
AN OPTIMIZATION PROBLEM 

A. Lagrangian shear for fronts 

We adapt the definition of Lagrangian shear in FBH— 
to the front propagation scenario. We first assume that 
u(r ,t) is given between an initial time to and a final 
time 1 1 = to + T, with no requirement that the flow 
be time-independent or periodic. Then the evolution of 
a front element over this time interval is obtained by 
solving Eqs. § We denote the flow map acting on the 
3D phase space by F = Fjf. Consider now an initial 
front R(s ) = (r(s),6(s)) at time to, and suppose the el¬ 
ement at a given position r(s) is displaced in the nor¬ 
mal direction an infinitesimal distance en. (The orienta¬ 
tion 8 of the front element is not perturbed.) Under the 



FIG. 2. (color online) The definition of Lagrangian shear p 
and normal repulsion p. An initial normal displacement en of 
a front element (left) evolves under F into a displacement ev 
whose tangential projection (blue) is —e times the Lagrangian 
shear p. The normal projection (blue) is e times the normal 
repulsion. 

time evolution F, the point R = (r, 8) maps forward to 
F{R) = Ri = (it,#!), and the relative displacement en 
between the original and perturbed front elements maps 
forward to a vector ev that will in general have a nonzero 
tangential projection onto the time-evolved front. (See 
Fig# We define the Lagrangian shear p over the time 
interval [to, U] to be this projected length divided by e. 

To compute p, we first note that the 3D normal dis¬ 
placement vector [n, 0] at to evolves forward to VF[n, 0] 
at ti, where VF is the 3x3 gradient of F, with com¬ 
ponents (VF)jj = F h j. Since v, shown in Fig. [2j is the 
ccp-component of VF[n, 0], we have 

v = 11., v V/ 'l Ij.yii, (3) 

where 

TT ( 1 0 0 \ 

Uxy ~ \ 0 1 0 ) 

is the restricted projection operator from the xyd-t&ngent 
space into the xy -tangent plane. Note [n, 0] = n^n. We 
then define 

[VF] xy = TL xy V FIll y ( 5 ) 

to be the 2x2 restriction of VF to the xy -tangent plane. 
Note also that Eq. ([I]) implies 

n = -fig = T^r'/lr'l, (6) 

where 

^((-o 1 )- p> 

Then Eq. (J3| becomes 

v = [X7F] xy n (8) 

= T[VF] xy Qr , /|r , |. (9) 

We next need an expression for the tangent vector gi to 
the front at iq. The 3D tangent vector to the curve R(s) 
at time to is of the form [g, a]. This vector maps forward 
under VF to a vector whose xp-component H xy VF(g, a) 
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is proportional to gi. Note that this fact is actually true 
regardless of the value of a, since a depends on the cur¬ 
vature dd/ds of r(s) and not its tangent. We may thus 
set a = 0 to obtain 


= [X7F\ xy g _ [VF] 

xy r 

\[VF} xy g\ |[VF] xy r'|’ 


where we again used the front-compatibility criterion 
Eq. |l]). We are now in a position to compute p via 


P = <-gi,v) = 


<[VF] xy r', [VF] xy Hr'} 


|r'||[VF] xy r' 


(r ',[VFff y [VF]^) 
^(r',r')(r',[VF]^[VF] xy r')' 


( 11 ) 


The minus sign on gi conforms to the convention of 

FBlM 

To summarize, we may express the Lagrangian 

shear as 


p{r,6) =p(r,r') 


(r\ C(r, 9)Q r') 
V( r ^ rl )( r ^C(r,9) r') 
(r',D(r,0) r') 
V(r',r')(r',C(r,0)rO’ 


where 


( 12 ) 

(13) 


C(r,0) = [VF£ y [VF] xy , (14) 

D(r,e) = ^(C(r,e)fl-nG(r,e)). (15) 

Note that C and D are both 2x2 symmetric matrices. We 
call C the projected Cauchy-Green tensor. Note that even 
though C is a 2 x 2 matrix acting on the rry-tangent plane, 
it still depends on the orientation coordinate 6. In the 
limit that the burning speed Vq vanishes, the projected 
Cauchy-Green tensor reduces to the usual Cauchy-Green 
tensor defined for passive advection. Finally, note that 
we have written p either as a function of (r, 9) or (r, r'), 
since p does not depend on the magnitude of r' and since 
9 is the angle of r', as given by Eq. (JTJ). 


B. Normal repulsion for fronts 

The bLCS will be defined to locally maximize the av¬ 
erage normal repulsion p. Adapting the definition of nor¬ 
mal repulsion in Refill to the front propagation scenario, 
we define p to be the normal component of the time- 
evolved displacement v in Fig. [2j That is, the normal 
repulsion is a local measure of the degree to which a 
nearby front is pushed away, with p > 1 indicating re¬ 
pulsion and p < 1 attraction. To compute p, we shall 
need the normal ni at the final time t\, which we claim 
is given by 


To prove this, we need only show that (gi, hi) = 0. From 
Eqs. (10 1 and (16), (gi,hi) is proportional to 

<[VF] xy r', [VF]- y 1T h) = ([VF]^[VF] xy r',h> 

= <r»=0. (17) 

Using Eq. ([ 8 ]), the normal repulsion is 


P= (ni,v) = 


([VF]- ir h, [VF] xy n) 
|[VF] xy 1T n| 


\j ([VF] xy 1T h, [VF] xy 1T n) 1 n 


(18) 


In summary, the normal repulsion is a function on the 
3D phase space given by 


p(r,9) = p{ r,r') 


1 

\/h • C' -1 (r,0)h 


(19) 


C. Variational problem 


The variational problem is now stated as follows. We 
seek a curve 7 , parameterized as r(s), that makes the 
average Lagrangian shear 


p [7] = - [ p(r,r')ds 

a Jo 


( 20 ) 


stationary under infinitesimal variations, where the aver¬ 
age is taken over the interval sq < s < Sq + a, choosing 
the initial point to have parameter so = 0. The end¬ 
points of 7 are assumed fixed under the variations. Note 
that p(r(s),r'(s)) does not depend on how the curve 7 
is parameterized, since the numerator and denominator 
in Eq. (13) are both second order in d/ds. However, the 


average in Eq. (20) does depend on the parameterization 


through the line element ds. Different choices of the pa¬ 
rameterization will change the weighting of the average 
along 7 . The stationarity condition of P[y] should be 
understood not only in terms of variations in 7 but also 
in terms of variations of the parameterization s along 
7 , with the interval of integration [ 0 , a] held constant. 
(In contrast, one could use the euclidean line element 
-y/(r', r')ds for the average, in which case ^[ 7 ] would be 
independent of the parameterization of 7 .) We shall find 
it useful to change the parameterization of the station¬ 
ary curves, so we next summarize the reparameteriza¬ 
tion of trajectories in the Hamiltonian and Lagrangian 
formalisms. 


D. Trajectory reparameterization for Hamiltonian systems 


NF} xy ^n 
| [VF] xy T n|" 


( 16 ) 


Suppose (R(s),P(s)) is a particular trajectory for an 
s-independent Hamiltonian H( r, p) and let E be the 
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value of the (conserved) Hamiltonian along the trajec¬ 
tory. Here s plays the role of time. Define a new trajec¬ 
tory parameterization r by 


dr/ds = n(r, p), 


( 21 ) 


for some phase space function k. and denote the reparam¬ 
eterized trajectory by (R(t),P(t)). Then (R(t),P(t)) 
is a trajectory of the Hamiltonian 


H e { r, p) = 


1 


«(r,p) 


{H{ r,p) - E), 


( 22 ) 


where the (conserved) value of He along the trajectory 
is E = 0. 

To prove this fact, we need only demonstrate X = 
{X, He}, where X(r, p) is an arbitrary phase space quan¬ 
tity, { , } denotes the Poisson bracket, and the dot de¬ 
notes differentiation with respect to r. To this end, 


{X ,H e } = {X 1 k~\H-E)} 


(23) 


= {X,k _1 }(R - E) + k _ 1 {X,R - E} (24) 


= K _ 1 {X, H} = K^X' = X, 
where the second to last equality utilizes 
X' = {X,H}, 


(25) 


(26) 


and the final equality uses Eq. (21 1 . 


E. Trajectory reparameterization for Lagrangian systems 

Suppose R(s) is a particular trajectory for the s- 
independent Lagrangian L(r, dr/ds). We again let E 
equal the value of the corresponding Hamiltonian along 
the trajectory. We also define a new parameterization r 
by Eq. (21), except that k = k( r,r / ) is expressed as a 
function of r' rather than p. Then the trajectory R(r), 
referred to the new parameterization, is a trajectory of 
the Lagrangian 

L E {r,r) = —(L(r,r) + E), (27) 

«(r, r) 

where L and h are L and k re-expressed as functions of 
. 1 , 


r = —r . 

K 


(28) 


To prove this, we recall the Legendre transforms 

H=(p,r')-L, (29) 

H e = {p,t)-L e , (30) 

where the momentum p is 

p = dL/dr' = dL E /dv. (31) 


Substituting Eqs. (29) and (|30|) into Eq. (22) yields 
Eq. <1271). 


F. Reparameterizing stationary curves of the average 
Lagrangian shear 


Defining the Lagrangian L = p, Eq. (13) yields 
L(r,r') = 


(r',D(r,e) r') 


V(r / ,r')(r',C(r, 0 ) r') 


(32) 


Re-expressing this in terms of the transformed derivative 
r, Eq. p 8 ), yields 


L{ r ,r) = 


(r,D r) 




(33) 


which has the same functional form as Eq. (32), since the 
r derivative occurs at the same order (twojin both the 
numerator and denominator. We now choose 


i(r,r) = 


1 


V(tr)(r,C: 


so that Eq. (27) becomes 


L E ( r, r) = (r, D r) + ^(r,r)(r,Cr) 


(34) 


(35) 


Equation (31) implies (p,r') = (r dL/dr 1 ), which is 


just the derivative of L in the direction of increasing ra¬ 
dius |r'| in the r'-space. Investigating Eq. (32), however, 


we see that the dependence of L on r' is solely through 
the angle 9 of r' and not |r'|. Thus, (p, r') = 0. Further¬ 
more, since r' = nr, (p,r) = 0 as well. Thus, Eqs. (29) 
and (30) imply 


H = L. 


H e = L 


E ■ 


(36) 


Since E = 0, the value of L E along the stationary trajec¬ 
tory R(t) also vanishes. These results were obtained by 
FBLpH through an alternative derivation. 

Theorem 1 Every stationary curve 7 of the average La¬ 
grangian shear Pfq], Eq. (20), is also a stationary curve 

of 


Pe[j] = — [ L E (r,r)d- 
T f Jo 


(37) 


= ^ J Q ^ r)) dr, (38) 

where E is the (conserved) value of the Lagrangian shear 
p along 7 . The value of L E along 7 is also conserved 
and equal to 0. Note that the parameterization R(r) of 7 
that makes P E stationary is typically different from the 
parameterization R(s) 0/7 that makes P stationary. 


Perfect shearless fronts 


Equations (32) and (36) show that E = 0 is equiv¬ 
alent to (r, Dr) = 0. Furthermore, one can show that 
any curve satisfying (r, Dr) = 0 is a stationary curve of 
Eq. (38). Thus, we have the following. 
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Theorem 2 A curve 7 satisfying (r, D r) = 0, or equiv¬ 
alently p = 0, is a stationary curve of 


Poll] 


^ J {r,Dr)dr. 


(39) 


We call such curves perfect shearless fronts (analogous to 
the perfect shearless barriers of FBH—), or just shearless 
fronts for short. They are the focus of the remainder of 
this article. 

Note that since D is symmetric and traceless [see 
Eq. (15)], it is like a pseudo-Riemannian metric tensor; 
“pseudo” because it has signature (—1,1). However, the 
resulting metric fails to be given purely by a bilinear 
form, since D itself depends on 0. Such a metric is called 
a (pseudo-)Finsler metric. Perfect shearless fronts are 
thus light-like, or null, geodesics of this pseudo-Finsler 
metric. 


H. Computing perfect shearless fronts 


In Sects. |III C-[ |IIIF[ and |IIIG[ fronts were viewed as 
curves in the iry-position space for the purpose of the 
variational problem. Referring now to the full xyO -phase 
space, a shearless front (r(r), 0(r)) satisfies both the front 
compatibility criterion Eq. (JT|) and (r, Dr) = 0. The lat¬ 
ter implies that (r, Cflr) = 0, and hence that r is an 
eigenvector of C( r, 0). One possible approach to finding 
the shearless fronts, then, is to integrate an eigenvector 
field of C(r,0). Let £(r, 0) denote a choice of unit eigen¬ 
vector over xyO- space. Taking r proportional to f. the 
front compatibility criterion Eq. 0 becomes 

xM) = 0, (40) 

where 

y(r,0) = n(0) • £(r, 0) = (sin0x — cos0y) • £(r,0). (41) 


This quantity is zero at (r, 0) if and only if g(0) = 
±£(r, 0). If the choice of £(r, 0) were smooth everywhere, 
then Eq. (40) would define a smooth two-dimensional 


constraint manifold in xyd- space. In general, however, 
there is no continuous choice of £(r, 0) available over 
the entire phase space. The continuity of £(r, 0) can 
break down at points in xyO -space where C( r, 0) has 
degenerate eigenvalues, and hence where there is not a 
unique eigendirection. Since degeneracy of real symmet¬ 
ric two-by-two matrices is a codimension-two criterion, 
the degeneracy points generically occur along curves in 
r0-space. Even if we remove these degenerate curves from 
consideration, £(r, 0) can not in general be chosen contin¬ 
uously on the remainder of .Ty0-space. This is due to the 
fact that £(r, 0) may obtain an overall rotation by 7t/2 
when transported along a loop encircling a degenerate 
curve. 


Rather than finding the eigenvectors explicitly, an im¬ 
proved approach is to use the constraint equation 


,0) = 0, 


where 


V’( r,0) = n(0) • C(r,0)g(0) 

= (sin 0x — cos 0y) • C(r , 0) (cos 0x 


(42) 


(43) 

■ sin 0y). 

(44) 


This quantity is zero at (r, 0) if and only if g(0) is 
any eigenvector of C(r,0), not just a particular choice 
£(r, 0) as above. The constraint surface Eq. (42) can 


thus be understood as the union of the two surfaces 


defined by Eq. (40) for the two choices of eigenvector 
f. The constraint surface Eq. (42) is a smooth sur¬ 


face without boundary (though potentially with self¬ 
intersections), which we call the shearless surface. 

A perfect shear less front, lifted into iri/0-space via the 
front compatibility condition, must lie within the shear¬ 
less surface. A perfect shearless front is thus an integral 
curve of a vector field tangent to the shearless surface. 
This vector field, normalized to unity, is 


g, 


where 


dr a 

d\ \/a 2 + b 2 
dd_ _ b 

dX \Ja 2 + b 2 ’ 


a = -^’ 
b = g • V0. 


(45) 


(46) 


An integral curve of (45) will lie within the shearless sur¬ 
face so long as its initial condition does. Note that these 
vector fields do not require us to find the eigenvectors of 
C(r,0). 

We shall also need the vector field that is everywhere 
orthogonal to Eq. (45), but still tangent to the shearless 
surface. This field is obtained by replacing g by n in 


Eqs. (45) and (46), i.e. 


dr a 

dn \J a 2 + c 2 n ’ 

dd c 


(47) 


dn yj a 2 + c 2 ’ 

where k parameterizes the curve and 

d , 

c = n • V«/j. 


(48) 


I. Burning Lagrangian coherent structures (bLCSs) as 
shearless fronts of maximal normal repulsion 


The normal repulsion p of a front element (r, 0) is given 
by Eq. (19). On the shearless surface, where n is an 
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eigenvector of C, this simplifies to 

p{ r, 0) = y/h- C(r, 9)h = vXl, (49) 


where Ax is the eigenvalue of C normal to the tangent. 
We now define the burning finite-time Lyapunov expo¬ 
nent (bFTLE) field as 

A = ^ log(Aj_) = ^ log(p). (50) 


Note that the specification of the eigenvalue A_l is based 
on the orientation of its eigenvector (normal to the front) 
rather than its magnitude; thus Ax need not be the 
largest eigenvalue. 

Since the bFTLE is defined on the two-dimensional 
shearless surface, we have overcome part of the chal¬ 
lenge of dimensions highlighted in Sect. |T| In particular, 
we can plot and visually interpret images of the bFTLE 
(Fig. |4^i). In some cases a bLCS follows a ridge of the 
bFTLE field. However, we see no reason why a ridge 
of the bFTLE field need satisfy the front-compatibility 
criterion, and so we opt not to compute bFTLE ridges. 
Again, we see the advantage of using the shearless fronts, 
for which the problem of defining the bLCS becomes 
the problem of selecting the “best” curve from the one- 
parameter family of shearless fronts. We choose the re¬ 
pelling bLCS to be that shearless front that (locally) 
maximizes the average normal repulsion 

R h] = \j C(r,6)ri dt, (51) 


where di is the xy-euclidean line element and L is the 
euclidean length of the curve. Similarly, the attracting 
bLCS is that shearless front that (locally) minimizes the 
average normal repulsion. 

There is still a question of what interval length L to 
choose. In keeping with the intuitive idea of following an 
bFTLE ridge, we shall adopt the approach of RefspSEH 
and restrict the length of a shearless front to those regions 
of the shearless surface where the second derivative of 
the bFTLE field in the direction normal to the front, 
given by Eqs. (47), is negative. In practice, the exact 
length used should make little difference, so we use the 
curvature criterion to tell us approximately where to cut 
the shear less curves. 


J. Cusps in the shearless fronts 

The shearless fronts, and hence a bLCS, may exhibit 
cusps when plotted in xy-space. This is one consequence 
of the 3D phase-space geometry that does not occur for 
traditional LCSs in passive advection. Cusps have al¬ 
ready been recognized as important features of BIMs for 
time-independent and time-periodic flows (see Fig. [3^,); 
cusps change the one-way barrier direction of the BIMs. 
Here, we determine when a cusp occurs along a shear¬ 
less front. First, we consider the fold of the shearless 


surface under the projection from the xy9 -space to the 
xy- space; a fold is where two branches of this projection 
meet. Since a shearless front is smooth on the shearless 
surface, if it remains on one branch of this surface, then 
its projection to xy-space is also smooth. Thus, a cusp of 
a shearless front can only occur at a fold of the shearless 
surface. What is less obvious, but nevertheless true, is 
that every time a shear less front passes around a fold, 
its xy-projection forms a cusp. To prove this, consider 
a shearless front (r(s),9(s)) as it passes around a fold. 
Locally, r can be expressed as a function of 9 , which by 
Eqs. (451 satisfies 


dr a „ 
de = b s ' 


(52) 


At the fold of the shearless surface, the 3D normal vector 
to the surface, equal to (V 0, dip / dO ), has no component 
in the 9 direction. That is, di^/d9 = 0. This implies that 
at a fold, a = 0 by Eq. ( |46| , and hence dr/d9 = 0, by 
Eq. (52). More precisely, one sees that, at a fold, dr/d9 
transitions from pointing in one direction to pointing in 
the opposite direction. This is exactly the criterion for 
forming a cusp—the curve r (9) moves forward and then 
backs up. In summary, then, a cusp occurs at exactly 
those points where a shearless front passes transversely 
through a fold in the shearless surface. 


IV. A TEST CASE 
A. Double-vortex in a windy channel 

As an example, consider an infinitely long channel 
flow containing two side-by-side, counter-rotating vor¬ 
tices plus a spatially uniform “wind” flowing from right 
to left. The streamfunction is 

i>{x,y) = ( 7 / 7 r)xexp(-x 2 )sin( 7 xy) + v w y , (53) 

where 7 is the vortex strength and v w is the wind speed 
(Fig. §• Throughout our analysis we use 7 = —1.0, 
v w = —0.15, and burning speed vq = 0.1. The channel 
width spans y = 0 to y = 1. We choose this flow because 
it models what can be realized in the laboratory of Tom 
Solomon®]. Furthermore, the BIM spans the channel 
without a cusp, which simplifies our analysis. Though 
in this paper we keep the wind speed constant, a future 
study will consider time-varying wind. Importantly, such 
time-varying wind can be realized in the Solomon lab 25 . 

For a steady wind, there are exactly two hyperbolic 
advective fixed points of the flow, assuming the wind is 
not too large. In Figs. and Sb these fixed points lie 
within the gray “semicircles” attached to the top and bot¬ 
tom channel wall, near the channel midpoint. The gray 
regions are “slow zones”, regions where the fluid speed u 
is less than the burning speed Vg. The bottom fixed point 
in these figures generates an unstable manifold transverse 









FIG. 3. (color online) Double vortex channel flow with 
vortex strength 7 = — 1.0 and burning speed Vo = 0 . 1 . a) no 
wind, v w = 0; b) v w = —0.15. The streamlines of the flow are 
illustrated by grey curves. The grey regions are “slow zones”, 
where the fluid speed u is less than vo- The blue and red 
curves are stable and unstable BIMs respectively. The BIMs 
begin at burning fixed points on the channel boundary. For 
plot a, the BIMs have cusps where they intersect the central 
slow zones. For plot b, the BIMs span the channel without 
forming a cusp. The arrows normal to the BIMs illustrate the 
burning direction. 


to the channel wall, while the top point generates a sta¬ 
ble manifold. The burning dynamics induces the bot¬ 
tom advective fixed point to split into two SSU burning 
fixed points, one at each intersection point between the 
boundary of the slow zone and the channel wall. (Only 
the rightmost burning fixed point is shown in Fig. !) 
There are also two SUU points near the bottom of the 
channel that do not concern us here. The right-burning 
SSU fixed point generates an unstable BIM, shown in red 
in Fig. [3j Similarly, the advective fixed point at the top 
of the channel splits into two SUU burning fixed points 
on the upper channel boundary. The leftmost of these is 
shown in Fig. [3j along with the attached stable BIM in 
blue. 

For weak or no wind, each BIM forms a cusp before 
crossing the channel (Fig. 0' .) As the wind speed in¬ 
creases above v w = Vq, each BIM crosses the channel 
without forming a cusp (Fig. s- .) The process by which 
the BIM attaches to the boundary is described in detail 
in RefP. For v w > vo, each BIM divides the channel 
in two. If the entire region left of the unstable BIM in 
Fig # is burned, then it will remain burne d for all time, 
forming a frozen front along the BI\ l 24 Fi l In contrast, 
for v w < Vq, any initially burned region will eventually 
propagate arbitrarily far down the channel to both the 
left and the right. Thus, the existence of an unstable 
BIM that traverses the channel with no cusp has a very 
clear experimental signature. 

A stable BIM traversing the channel without a cusp 
(Fig. I^)) forms a kind of basin boundary. An initial 
point stimulation left of the stable BIM will eventually 
grow to become a frozen front. However, a stimulation 



FIG. 4. (color online) a) The two-dimensional shearless 
(ip = 0 ) surface embedded in the three-dimensional phase 
space. The coloring represents the two-dimensional bFTLE 
field on this surface. The Cauchy-Green integration time is 
T = 3. The black curve is the stable BIM attached to the 
rightward propagating burning fixed point at the top of the 
channel. The computation of the isosurface is restricted to 
the computation box shown, b) A view of the bFTLE field 
on the shearless surface projected onto the rry-plane. The 
two white patches on the bottom left and right are where the 
shearless surface has exited the computation box. c) The sign 
of t he b FTLE curvature normal to the shearless fronts (see 
Eq. (471. Blue is negative and red is positive. 


right of the stable BIM will eventually fill up the entire 
channel. In the remainder, we shall focus on the stable 

BIM. 


B. The bLCS for a steady wind—short integration time 

We now apply the method of Sect. [TTT] to the steady 
wind case (v w = —0.15) with the relatively short Cauchy- 
Green integration time T = 3. We should find a bLCS 
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that follows the initial length of the time-independent 
BIM. Before finding the bLCS, however, it is helpful to 
visualize the shearless surface. Fig. shows the surface 
in 3D. For comparison, the stable BIM is shown in black. 
Near the burning fixed point, the BIM lies near (though 
not precisely within) the shearless surface. This is an 
important first check on the validity of our technique; if 
the BIM were far from the shearless surface, we would 
have no chance of recovering it. 

At the bottom center of Fig. |4jr is a “pitchfork” in¬ 
tersection between two branches of the surface. The dif¬ 
ferentiability of the if) field implies that the intersection 
curve between these two branches must be a shearless 
front, i.e. a trajectory of Eqs. (45). Thus, no shearless 
fronts may pass through this intersection. 

The coloring in Figs. [4^, and[4]D represents the bFTLE 
field Eq. (50). Notice that the BIM follows a bFTLE 
ridge at the top of Fig. [4j3, another critical check on our 
technique. Curiously, there is a second bFTLE ridge on 
the left in Fig. This ridge lies on a separate branch 
of the pitchfork structure, as seen in Fig. [4ji. 

We next identify those regions where the second deriva¬ 
tive of the bFTLE field normal to the shearless curves is 
negative. Figure [4}' illustrates the sign of this curvature, 
with blue representing negative and red representing pos¬ 
itive. There are two regions of negative curvature, one 
for each of the bFTLE ridges. We focus on the rightmost 
region associated with the BIM. 

We next compute shearless fronts, targeting the bF¬ 
TLE ridge on the right. We select initial conditions using 
the intersection of the shearless surface with a plane nor¬ 
mal to the BIM, yielding the bold line of initial points in 
Fig. [5^P^. We integrate Eqs. |45j upward toward the top 
of the channel and downward toward the bottom. Fig¬ 
ures and[5}:> show these shear less fronts in 3D and 2D, 
respectively. 

We then compute the average normal repulsion, 
Eq. ( [5l] ) , along the segments of the shearless curves above 
the dashed line y = 0.7 in Fig. [5]:). This line is chosen to 
lie close to the bottom of the negative curvature region in 
Fig- U- The resulting average normal repulsion is shown 
in the inset of Fig. [5jc. It exhibits a clear local maximum, 
which selects the red shear less curve isolated in Fig. §5- 
This curve is the bLCS, and it faithfully follows the ini¬ 
tial length of the BIM, as anticipated, before eventually 
deviating markedly from it. 




-0.1 0 0.1 0.2 


FIG. 5. (color online) a) Shearless fronts computed for 
Cauchy-Green integration time T = 3. The initial points 
(bold) of the trajectories were computed by intersecting a 
plane transverse to the BIM (at the red dot) with the shear¬ 
less surface. The BIM is the thin black curve, b) The shearless 
fronts are shown projected onto the xy- plane, c-inset) Plot of 
the average normal repulsion, Eq. (511, as a function of the 
x-coordinate of the initial point for each trajectory in parts a 
and b. A smooth maximum is marked in red and determines 
which shearless front is the bLCS. c) The bLCS is plotted 
in red next to the BIM in black. The open circle marks the 
initial condition for generating the bLCS. 


V. CONCLUSION 


C. The bLCS for a steady wind—long integration time 

We extract a bLCS for the same dynamics as in 
Sect. |IVB| but for a longer Cauchy-Green integration 
time of T = 6. The result is summarized in Fig. [6] 
Clearly the bLCS (red) closely tracks the BIM (black), 
considerably longer than for the T = 3 case. In fact, the 
agreement is quite good all the way to the bottom of the 
channel. This is exactly what we expect to see as the 
integration time T is increased. 


We have developed a technique to extract maxi¬ 
mally repelling (or attracting) fronts, called burning La- 
grangian coherent structures (bLCSs) for the finite-time 
dynamics of fronts evolving within a flowing fluid. We 
verified that this technique returns the burning invariant 
manifolds (BIMs) for the case of a time-independent flow. 
Clearly more work is needed to establish the relevance of 
these bLCSs to time-periodic and, most importantly, to 
time-aperiodic flows. This will be pursued in a separate 
publication. 
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FIG. 6. (color online) The bLCS (red) and the BIM (black) 
for Cauchy-Green integration time T = 6. The open circle on 
the bLCS marks the initial point used to integrate the shear¬ 
less front. (The neighboring shear less fronts are not shown.) 
The inset shows the average repulsion on a shearless front as 
a function of the ^-coordinate of its initial point. The average 
was computed along the shearless fronts above the dashed line 
at y — 0.2. 
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